Global stability and bifurcations in a mathematical model for the waste plastic management in the ocean

The use of plastic is very widespread in the world and the spread of plastic waste has also reached the oceans. Observing marine debris is a serious threat to the management system of this pollution. Because it takes years to recycle the current wastes, while their amount increases every day. The importance of mathematical models for plastic waste management is that it provides a framework for understanding the dynamics of this waste in the ocean and helps to identify effective strategies for its management. A mathematical model consisting of three compartments plastic waste, marine debris, and recycle is studied in the form of a system of ordinary differential equations. After describing the formulation of the model, some properties of the model are given. Then the equilibria of the model and the basic reproduction number are obtained by the next generation matrix method. In addition, the global stability of the model are proved at the equilibria. The bifurcations of the model and sensitivity analysis are also used for better understanding of the dynamics of the model. Finally, the numerical simulations of discussed models are given and the model is examined in several aspects. It is proven that the solutions of the system are positive if initial values are positive. It is shown that there are two equilibria \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E^0$$\end{document}E0 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E^*$$\end{document}E∗ and if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal {B}}}{{\mathcal {R}}}<1$$\end{document}BR<1, it is proven that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E^0$$\end{document}E0 is globally stable, while when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal {B}}}{{\mathcal {R}}}>1$$\end{document}BR>1, the equilibrium \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E^*$$\end{document}E∗ exists and it is globally stable. Also, at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal {B}}}{{\mathcal {R}}}=1$$\end{document}BR=1 the model exhibits a forward bifurcation. The sensitivity analysis of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal {B}}}{{\mathcal {R}}}$$\end{document}BR concludes that the rates of waste to marine, new waste, and the recycle rate have most effect on the amount of marine debris.

parameter µ is the recycled waste rate to be reproduced as new waste and finally by θ we denote the recycled waste rate to be lost.
Using the properties of the Levenberg-Marquardt backpropagation (LBMBP), authors in 1 designed an artificial neural network method for solving this model computationally.In 2 , the model was solved using two schemes based on a modification of the Morgan-Voyce (MV) functions by applying directly matrix collocation procedure and by using the quasi-linearization together with the modified MV collocation method.As mentioned, in previous studies, this model has been investigated numerically, but so far, its dynamics have not been comprehensively discussed in any work.In the present study, we intend to analyze this model with mathematical tools, although we also test the obtained results by simulation and numerical examples.
Compartmental models have been used in mathematical modeling to study the epidemics and ecological systems [3][4][5][6] .Also, there are some three-compartmental epidemic models to study diseases, including SIR or SIS models with vaccination or quarantine [7][8][9][10][11] , which may be useful for studying this model due to the nature of the cycle of propagation and recycling of plastic waste in the ocean.But according to our knowledge, the system studied in none of the works are completely compatible with the system related to this phenomenon in this paper.
In the present work, we first give some basic properties of the model in Section "The description of the model".After showing that the solutions of the system are always non-negative, we obtain two equilibrium points for the system along with their existence conditions.Then, by applying the next generation method, we determine the basic reproduction number BR for the model.The investigation of the long-term behavior of the model as a dynamical system is done through the analysis of the WPM system (1) stability at the equilibrium points in Section "Stability of the model".Also, we state and prove the conditions in terms of this quantity under which the system is stable at each equilibrium point.The bifurcations of the model and the sensitivity analysis of the parameters are discussed in Sections "Bifurcations of the model" and "Sensitivity analysis", respectively.Numerical experiences and simulation of the solutions of the system are also carried out via some examples in Section "Numerical experiments".Finally, the results of the paper are summarized as conclusions.

The description of the model
In this section we obtain some basic properties of the model such as the invariant set of solutions, the equilibria of the model and the basic reproduction number.

Non-negativity of solutions
The following lemma states that the solutions of system (1) remain non-negative with positive initial values.
Proof Assume W(0) > 0 , M(0) > 0 and R(0) > 0 .We let and so τ > 0 since all variables are continuously differentiable.Now, if τ = +∞ then the positivity of solutions holds, but if 0 < τ < +∞ then one of the variables is zero at τ and for t > τ is negative.Assume that W(τ ) = 0 and W(t) < 0 for t > τ .From the differential equation corresponding to W we find Therefore W is non-decreasing at τ and W(t) ≥ 0 for t > τ .This is a contradiction with the previous assump- tion for τ < +∞ and thus for all t we have W(t) ≥ 0 .We can use a similar argument for variables M and R and conclude M(0) ≥ 0 and R(0) ≥ 0 for all t > 0 .Therefore the solutions of the system (1) are always non-negative and � = {(W, M, R) : W ≥ 0, M ≥ 0, R ≥ 0} is positive invariant for solutions of system (1).

Equilibria of the model
Solving the following system of equations gives the equilibrium points (in the format of ( M, W, R) ) for model (1): Two equilibrium points are obtained; the equilibrium if M = 0 , and the equilibrium when M > 0. τ = sup{t > 0 : W(t) ≥ 0, M(t) ≥ 0, and R(t) ≥ 0}, Vol Therefore if we let BR = β (µ+θ) αγ θ , then M * = γ β (BR − 1) > 0 if and only if BR > 1 and thus E * exists if and only if BR > 1 .On the other hand, considering that the components of E 0 are always non-negative, we have stated the following lemma.
Lemma 2 For model (1), the equilibrium point E 0 always exists, while the equilibrium point E * exists under the condition that BR > 1.
We call BR as the basic reproduction number of the WPM model.The basic reproduction number for a compartmental model is actually defined as the number of secondary reproductions in a completely homogeneous population by a member of a compartment that causes the spread of undesirable properties in other compartments 12 .In the preceding discussions, this quantity was implicitly given by using the condition that the amount of marine debris ( M * ) is positive at the equilibrium point E * .Moreover, it can be also obtained by the next generation matrix method as the spectral radius of the Jacobian matrix of the model at E 0 as follows 13 .
In this method we rewrite the second equation in (1) corresponding to the marine debris (compartment M) as where F = βMW and V = αM .Here, F refers to the terms in the equation that generate new amounts of marine debris (M), and V consists of terms that represent the transmissions between compartment M and other compartments.

Stability of the model
In the paper 2 , the authors studied the local asymptotic stability of the model.By examining the eigenvalues of the Jacobian matrix at the equilibrium points of the model and using the Routh-Hurwitz criterion 14 , they determined the conditions under which the eigenvalues of the Jacobian matrix have a negative real part and showed that when BR < 1 , the system is stable at equilibrium point E 0 and if BR > 1 , this point is unstable and equilibrium point E * is stable.Now, in the following we investigate the global stability of the model at the equilibria of the system (1).

Theorem 1
The equilibrium E 0 of model ( 1) is globally asymptotically stable when BR ≤ 1.
Proof For this purpose, we use the basic concept of global stability that states under the given conditions, all solutions of the model converge to the equilibrium point E 0 starting from any point 15 .
From third equation in system (1) we have The first equation in (1) implies W ′ (t) ≤ − γ W(t) + µR(t) and thus we get = W 0 and lim t→∞ R(t) = θ = R 0 .Therefore all trajectories of (4) converge to (W 0 , R 0 ) and the equilibrium E 0 is globally asymptotically stable.But, if µ + γ < θ we have The next theorem states the conditions under which the equilibrium E * is globally stable.For this purpose we use the Lyapunov's direct method 16 which has also been employed by many authors [17][18][19][20][21] Theorem 2 The equilibrium E * of model ( 1) is globally asymptotically stable if BR > 1 and α = γ.
Proof Using combinations of composite quadratic and common quadratic terms 17 , we consider function Function L is C 1 on the interior of defined domain set, and we also see L ≥ 0 , the equilibrium point E * is the global minimum of L, and moreover L(E * ) = 0 .By differentiating L we get From system (2) we have Thus we can write and thus since γ = α , we have Also, dL dt = 0 if and only if R = R * (which implies from (1) M = M * and W = W * ) and Therefore the equilibrium E * is globally asymptotically stable according to LaSalle's invariant principle 16 .

Bifurcations of the model
We now consider the behavior of the model when BR = 1 .For this purpose we choose β as the bifurcation parameter and from BR = 1 we find To identify the bifurcations that model (1) exhibits, we apply center manifold theory 22 .According to this method we need to carry out the following change of model variables.Let x 1 = W , x 2 = M and x 3 = R by using vector X = (x 1 , x 2 , x 3 ) T , the waste plastic management model can be rewritten in the form X ′ = F(X) The jacobian matrix is obtained as x and at equilibrium E 0 for β * (when BR = 1 ) is because β * (µ+θ) γ θ = α .The eigenvalues of J 0 β * are 1 = 0 and those for the sub-matrix For this matrix we obtain and thus its eigenvalues have negative real part 14 .
The right eigenvector u = (u 1 , u 2 , u 3 ) T of matrix J 0 β * corresponding to the zero eigenvalue is obtained by solving the system J 0 β * u = 0 and we get u 1 = −α , u 2 = γ and u 3 = 0 .Besides, the components of the left eigenvector v = (v 1 , v 2 , v 3 ) corresponding to the zero eigenvalue can be found from v u = 1 as are gotten as v 1 = 0 and v 2 = 1 γ ad v 3 = 0 .It can easily be checked that vJ 0 β * = 0 .According to the center manifold theory we have to calculate the two following coefficients All second order partial derivatives of f k , (k = 1, 2, 3) with respect to x i , (i = 1, 2, 3) and β are calculated at γ θ , 0, θ .Since v 1 , v 3 and u 3 are zero, we only need to calculate the fol- lowing second order partial derivatives: Thus Therefore, according to Theorem 4.1 in 22 since a < 0 and b > 0 we conclude that there is a forward bifurcation (transcritical bifurcation) at β * (when BR = 1 ) for the waste plastic management model.

Sensitivity analysis
To find out how sensitive a function is to changes in the variables in its formula, the method of sensitivity analysis is used.This method uses a quantity called the normalized forward sensitivity index as a measure of the sensitivity of each variable.Since the basic reproduction number has a significant impact on the behavior of the model, we calculate the sensitivity indices of BR for its variables, which are β , , µ , α , γ , and θ.
The normalized forward sensitivity index of variable BR for a variable v, is defined by If SI BR v > 0 then the variable v has positive impact on BR and the value of variable BR increases by increasing the value of v.While SI BR v < 0 shows the reverse impact of v on BR ; increasing the value of v implies decreasing in the value of R 0 .Also, the magnitude of SI BR v shows the proportion of changes in BR with respect to v.
The normalized forward sensitivity indices for BR are caculating as follows: We see that SI BR v > 0 for v = β, and µ , but SI BR v < 0 for v = α, γ and θ .Thus, any increase (or decrease) in values of β, and µ has direct impact on the value of R 0 .However, any increase (or decrease) in values of v = α, γ and θ has reverse impact on the value of BR.

Numerical experiments
We consider the values β = 0.15, γ = 0.41, α = 0.65, µ = 0.4, = 0.36, θ = 0.15 for the parameters in model (1).For these values we have BR = 0.743 < 1 and according to Theorem 1 the equilibrium E 0 = (0, 3.2195, 2.4000) is stable.Figure 1 illustrates the solutions of the model with these parameter values and 20 different initial values for sub-populations.Moreover, the left picture in Fig. 3 shows the solutions of the system with initial values as W 0 = 1.5 , M 0 = 2 , and R 0 = 1.Now, we change the parameter values to β = 0.4, γ = 0.21, α = 0.5, µ = 0.4, = 0.66, θ = 0.2 .Here, we have BR = 7.5429 > 1 and by Theorem 2 the equilibrium E * = (1.25,3.44, 3.30) is stable.In Fig. 2 the solutions of the system (1) has been depicted for 20 initial values for sub-populations.In the case of that the initial values are supposed as W 0 = 1.5 , M 0 = 2 , and R 0 = 1 , the right picture in Fig. 3 shows the solutions of the model.
For parameter values γ = 0.41, α = 0.65, µ = 0.4, = 0.36, θ = 0.15 (as in the first example), if we solve the equation BR = 1 for parameter β as the bifurcation parameter, we get β * = 0.2019 .Thus for this value the model exhibit a forward bifurcation as it can be seen in bifurcation diagram presented in Fig. 4. Now, we investigate the impact of the parameters of the model on the dynamics of the waste management system by using sensitivity analysis of the basic reproduction number BR with respect to each parameter as it was explained in Section "Sensitivity analysis".Let us consider the parameter values in Table 1 for parameters in the model introduced by system (1) and their corresponding sensitivity indices with respect to BR as a differentiable function.
The normalized sensitivity indices for parameters have been shown as a chart in Fig. 5.According to the Table 1 we find that the sensitivity indices for parameters β , and µ have positive value and their impact on BR is direct.While, the parameters α , γ and θ have negative sensitivity indices and thus they have reverse impact on BR .Thus for example, an increase (decrease) in values of and µ by %10 yields to a %10 and %8.889 increase (decrease) in BR , respectively.For instance, for parameter values in Table 1 we have BR = 42.4286 and if we decrease by %10, the basic reproduction number also decreases by %10 and becomes BR = 38.1857 .On the other hands, a %10 increase (decrease) for example in γ and θ causes a %10 and %8.889 decrease (increase) in BR .For example, if the value of θ is increased by %20 (to θ = 0.06 ), then the value of BR decreases by %17.778 and becomes BR = 34.8763 .Therefore, according to Table 1 we find that either decreasing the rates of waste to marine ( β ) and new waste ( ) or increasing the recycle rate ( α and γ ), have the most impact in reducing the value of BR and as a result most impact on controlling the amount of marine debris.The impact of parameters β and α on the model have also been depicted in Figs. 6 and 7, respectively.The parameter values are assumed as β = 0.75, γ = 0.21, α = 0.5, µ = 0.4, = 0.66, θ = 0.05 and initial values are W 0 = 1.5 , M 0 = 2 , and R 0 = 1 .In Fig. 6 the values of β change in interval [0.1, 2.1] and it is seen that by increasing β ,the final solution cor- responding to the amount of the marine debris (M) leaves the zero and will take positive values.Indeed, the stability of the model changes from E 0 to E * .With the same terms and for values of α in [0.1, 0.9], the solutions of the model have been shown in Fig. 7.By increasing the value of α , we observe that the values of M finally take zero value.This shows that by increasing α the equilibrium point E * becomes unstable and E 0 becomes E * stable.

Conclusions
In this paper, we studied the waste plastic management (WPM) system in the ocean through a mathematical three-compartmental model.The basic reproduction number BR , and two equilibria of the model were found, in addition to positivity of solutions of model.The dynamics of the model was determined in terms of threshold BR ; if BR < 1 , it was proved that the equilibrium E 0 is globally stable, while the equilibrium E * exists and it is stable when BR > 1 .Also, it was shown that the model exhibit a forward (transcritical) bifurcation at BR = 1 .
The sensitivity of the model has been analyzed by calculating normalized forward sensitivity index for each parameter for BR and it was concluded that decreasing the rates of waste to marine ( β ) and new waste ( ) or increasing the recycle rate ( α and γ ), are most effective for controlling the amount of marine debris.Finally, the theoretical results were discussed also numerically for different parameter values and various initial values for sub-populations via several examples, solutions of model and bifurcation diagram.The global stability of the equilibrium point E * has been proved for the case that the marine debris recycling rate ( α ) and the direct recycling rate of waste materials ( γ ) are equal.Constructing a more appropriate Lyapu- nov function that does not impose such an additional condition on stability can be the subject of future studies.Investigating the impact of seasonality on the system behavior may also complement the present study, since the waste rate to become marine ( β ) can be considered as a periodic function.

Fig. 3 .
Fig.3.Solutions of the model for two cases BR < 1 (the left picture) and BR > 1 (the right picture) with their theoretical solutions when initial values are same.
Basic reproduction numberFirst, we note that M * in equilibrium point E * can be written as Thus, the equilibrium point E * exists if and only if